Polyelectrolyte chains in poor solvent. A variational description of necklace 

formation. 

Gabriele Migliorini 1 Vakhtang Rostiashvili 1 Namkyung Lee 1 and Thomas A. Vilgis 1,2 
1 Max Planck Institute for Polymer Research 
10 Ackermannweg, 55128 Mainz, Germany. 
2 Laboratoire Europeen Associe, Institute Charles Sadron 
6 rue Boussingault, 67083 Strasbourg Cedex, France. 
(February 1, 2008) 

We study the properties of polyelectrolyte chains under different solvent conditions, 
using a variational technique. The free energy and the conformational properties of a 
polyelectrolyte chain are studied minimizing the free energy Fn, depending on N(N — 
l)/2 trial probabilities that characterize the conformation of the chain. The Gaussian 
approximation is considered for a ring of length 2 4 < N < 2 16 and for an open chain of 
length 2 4 < TV < 2 9 in poor and theta solvent conditions, including a Coulomb repulsion 
between the monomers. In theta solvent conditions the blob size is measured and found 
in agreement with scaling theory, including charge depletion effects, expected for the case 
of an open chain. In poor solvent conditions, a globule instability, driven by electrostatic 
repulsion, is observed. We notice also inhomogeneous behavior of the monomer-monomer 
correlation function, reminiscence of necklace formation in poor solvent polyelectrolyte 
solutions. A global phase diagram in terms of solvent quality and inverse Bjerrum length 
is presented. 



INTRODUCTION 



Polyelectrolytes are macromolecules with ionizable groups that can dissociate in solution to produce 
polyions (polymers with fixed or mobile charges) that arc cither positively or negatively charged and 
counter-ions of opposite charge in solution. Polyelectrolytes are found in many biological systems (pro- 
teins, nucleic acids, ..) as well as in synthetic products (resins, ..). As an example of polyelectrolyte 
we mention the polystyrene sulfonate {C^HjSO^,Na) n that dissociates in water to produce polyanions 
(C%H'jSO^)n and counter-ions Na + . The behavior of polyelectrolytes is much less understood systems 
compared to others commonly discussed in macromolecular science. The main reason is the difficulty in 
applying renormalization group and scaling ideas in systems where long ranged forces are present. More- 
over, the special features of the Coulomb potential below (long ranged) and above four (short ranged) 
dimensions prevent from a clear physical picture even in three dimensions. For weakly charged and 
flexible polyelectrolytes the situation is less severe First, the flexibility allows to use the Wiener 

measure, i.e., Gaussian chains in absence of interactions, and electrostatic corrections to the stiffness of 
the chain Q can be safely ignored. Second, in three dimensions the Coulomb potential is long ranged 
(I/7') and variational techniques are able to provide physically sensible results. 

Variational methods have been discussed in the context of polymer physics |^-^] and these statements 
have been carefully investigated, with the result that variational methods are well suited to study systems 
with long ranged interactions. In the last two decades considerable effort have been devoted to the 
statistical properties of polyelectrolyte solutions, 10-12 . It is usual custom to study strongly charged 
macromolecules, where the Coulomb interaction between the monomers is larger than the short range 
interactions of the van-der-Waals type. Here we will investigate both regimes by means of an extended 
variational method, including the case of weakly charged polyelectrolyte chains ll&O]. We are going to 
study the collapse of a single chain, when the short ranged interactions dominate and recover the limit 
of strongly charged chains, where elongated conformations are expected [ pd] , |l5[ . 

The main challenge is to study the case of poor solvent conditions and the formation of globular states. 
So far, these regimes are mainly studied by scaling considerations and blob pictures |i"6|-|lS|] as well as 
Monte Carlo simulations jl9| and by mean field approaches The essential features of a neutral 

chain in poor solvent can be summarized as follows. Attractive and repulsive interactions (two and three 
body terms in the virial expansion) balance each other; as a result a globule of size R ~ 6(/V/r) 1 / 3 for ms, 
where r defines the relative temperature distance from the theta point, r = T — 0|/0, and determines the 
second virial coefficient v — — b 3 r. The globule itself can be viewed as a melt of blobs of diameter £ = b/r. 
The density of the globule is then given by c = r/6 3 . The main difficulty, however, is to obtain these 
properties from a more rigorous point of view. This seems to be nontrivial, since most of the properties 
of the globule are determined by the surface tension 70 ~ k^T/^ 2 . However, variational techniques are 
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suitable for these cases, since fluctuations inside globules are small, and play only a significant role on 
length scales smaller than the thermal blob size £, where entropy dominate and the chain conformation 
is Gaussian. The last remark applies also for the surface of the globule. 

In the simple case of neutral chains, we investigate the coil to globule transition driven by the short 
ranged interactions. Globule formation, as a balance between the two and three body terms, takes 
place under poor solvent conditions. The physics of uncharged macromolecules has been nowadays well 
understood |ll],|7],^0|,|lj . Charged macromolecules are far less understood systems because of the interplay 
of short ranged attractive forces and long ranged Coulomb forces. It is our first aim to recover the scaling 
results for a neutral chain within the two variational schemes introduced in section II. Previous attempts 
to describe poor solvent chains via a variational method [p2| and the description of the necklace phase E| 
should be mentioned as well. Once the correct scaling behavior is recovered, we will return to the problem 
of a weak Coulombic perturbation and its effect on the macromolecule conformations. This is discussed in 
detail in section III and IV. As expected from standard energy considerations p5|-p7t the globule is stable 
for charges smaller than a critical threshold (in a similar fashion as a charged liquid drop). For charges 
exceeding this point, the globule will elongate to minimize the Coulomb repulsion. Due to connectivity 
of the macromolecule, the globule will split into two which are connected by a stretched string of thermal 
blobs. In the following we will explore this scenario within the variational method for an open chain in 
section IV. The properties of the necklace phase are discussed and a global phase diagram in the space 
of solvent quality and inverse Bjerrum length is presented. 



I. POLYELECTROLYTE PHENOMENOLOGY 



A. Polyelectrolyte Chain in Theta Solvent 



In this section we review the fundamental scaling results for polyelectrolyte chains in the case of a dilute 
solution below the overlap concentration, such that only single chain statistics can be taken into account 
jl5],^8|]. We assume first to be in the weakly charged limit, where the charge fraction is small enough to 
assume Gaussian statistics for the non interacting chain and to ignore counter-ion condensation. At high 
temperature, the interaction between the counterions and polyelectrolyte chains in solution is weak and 
the distribution of counter-ions is homogeneous. Only a negligible small fraction of counter-ions are in 
the neighborhood of the chain, and their influence on the chain conformation can be neglected. Chain 
entropy and Coulomb interaction will then determine a globally stretched conformation of the chain. 
At smaller scales, entropy dominates the behavior and the usual blob picture holds ^8j. The chain is 
made of a string of electrostatic blobs of size £ e i, each made of g \ monomers. The size of each blob 
and the number of monomers within one blob can be determined by simple energy considerations. The 
electrostatic energy of a blob is of the order of the thermal energy, i.e. (/3oi) 2 'b/£c1 — 1 and the Bjerrum 
length is defined as 1-q — e 2 /ek^T, where e is the dielectric constant of the solvent and T the temperature. 
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Moreover the polymer segments inside the blob obey Gaussian statistics, i.e. £ c i = bgj , so that the size 
of the blob is given by £ c i ~ 6(6//b) 1 ^ 3 / -2 ^ 3 - The length of the polyelectrolyte is then R ~ (N/g e \)£ e \, so 
that 

R ~ b(l B f/b)V 3 N (1.1) 

A numerical solution of the variational method described in section II, confirms the scaling behavior 
predicted by equation (|l . l|) , including logarithmic corrections in the form N(\ogN) 7 Moreover, we were 
able to derive equation (|1.1| ) by assuming the proper ansatz for the monomer-monomer correlation 
function at large distances, assuming a non universal dependence of the exponent 7 at small values of 
the interaction strength. Differently, a mere implementation of des Cloizeaux asymptotic analysis for 
polyelectrolytes would lead to an unphysical overstretching, i.e., R cx N 2 for the end to end distance. 



Details are described elsewhere [p|]29| ■ Presence of logarithmic corrections to equation (IT) were already 
suggested in the past |3(|, with a universal exponent 7 = 1/3. The presence of logarithmic corrections, 
that makes the blob size a non-universal quantity, appears naturally in the variational approach and finds 
its origin in the mean field nature of the method, similarly to logarithmic corrections in a field theory 
at the upper and lower critical dimension. The electrostatic blobs close to the end of the chain become 
larger, since end monomers experience less electrostatic interactions compared to monomers in the middle 
of the chain, leading to inhomogeneous fluctuations in the electrostatic potential. These effects have been 
recently studied by an analysis of the classical path of the polyelectrolyte chain |n|. The resulting 
generalized expression for the blob size reads 



2 



t = ^ n 2) 

[M^g^ + i]^ ' ( ' j 

where £ G i is the usual electrostatic blob size, obtained from scaling arguments. The quantities L and £m 
are respectively the length of the chain and the maximum blob size at the chain ends. It is important to 



remark that, neglecting the z dependence in equation (1.5 ) one recovers the variational result of de Gennes 
et al. pd| , with the proper logarithmic behavior. Shortly below, we show how these correlations can be 
detected as well, by a more subtle variational technique, which uses heterogeneous monomer-monomer 
correlation functions as trial function for the chain statistics. 



B. Polyelectrolyte Chain in Poor solvent 



We consider first a polyelectrolyte chain in poor solvent, where the monomer statistics becomes much 
richer, compared to that related to the theta solvent. The polyelectrolyte chains consists of N monomers 
of size b and the fraction of charged monomers /. The uncharged chain in a poor solvent forms a 
globule. The monomer density p — rb~ 3 inside the globule is defined by the balance between the two 
body attraction and the three body repulsion. From standard scaling arguments we can distinguish four 
different regions. At weak enough values of the solvent quality the chain is Gaussian and R — bN 1 / 2 . 
Below the critical value of r c = A -1 / 2 we enter the globular phase where the overall size of the chain 
scales as 

/jV\ 1/3 

R = b\-j ■ (1.3) 
We review briefly the scaling theory when Coulomb interactions are present to the poor solvent chain 



[|13| , [25 32 33j] . For values of the two body attraction lower than the critical threshold r < t c the chain will 
extend to the electrostatic blob chain discussed above if the electrostatic interaction exceeds IbJ 2 /b > 
N~ 3 / 2 , and therefore overcomes the Gaussian entropy. In the poor solvent regime r > t c the globules 
remain spherical for small charge fractions. However, if the charge fraction increases and exceeds the 
value 
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the globule changes shape and an elongated structure emerges, still maintaining the local scaling prop- 
erties of a dense globule. The extended globule is referred in the literature as a cylindrical globule J25|. 
In absence of charges the surface tension 70 of the globule is of the order of kT per blob at the globular 
surface, i.e. 70 « kTt;~ 2 . When the charge fraction becomes larger, the Coulomb repulsion between the 
monomers will finally dominate over the surface free energy therefore affecting the shape of the glob- 
ule. Whenever the Coulomb repulsion F cou \ k, k^Tl-BpN 2 / R becomes comparable to the surface energy 
F sur w 70-R 2 oc kBTN 2 / 3 r 4 / 3 , the critical Bjerrum length at which this occur can be simply derived by 
comparing the bulk and surface free energy of the globule given above p7[ . In this regime the average 
length of the cylindrical globule is given by 
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R\\ = b ~ ( • (1-5) 
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At higher values of the critical Bjerrum length defined by eq.( |l.4| ) the globule eventually split into two 
and more connected globules, according to a Rayleigh instability of charged spheres p6| , p7pq ]. At even 
larger charge fractions the chain finally stretches and the standard electrostatic blob regime is recovered, 
as soon as 

;** r2 
** l B J 



> r\ (1.6) 

the chain becomes extended and the electrostatic and thermal blob size become equal. The chain has 
then an extension R = bNr. The foregoing scenario is depicted in Fig.l for a ring and an open chain, 
neglecting for the moment a necklace formation, as an alternative, energetically favorable conformation to 
the cylindrical globule introduced by Khokhlov. Below, we investigate these different regimes predicted 
by scaling theory using a new variational technique we recently introduced S . It is important to remark 
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that connectivity plays here a crucial rule. All sort of scaling arguments do not take connectivity into 
account, while we will be able to describe the globule to necklace transition by mean of a variational 
method that starts from the original Hamiltonian for the problem and properly include the connectivity 
of the chain. We consider a generalized form of the Gaussian variational method j^J^] to study the 
scaling properties and spatial conformation of a polyelectrolyte globule in poor solvent, for different 
values of the Coulomb interaction strength between the monomers, obtaining an independent test for 
the results commonly discussed in the scaling theory approach. By mean of the generalized Gaussian 
variational method we will be able to study the overall extension of the globular structure measuring, 
e.g. the end to end distance and the internal properties of the extended globule, e.g. point to point 
correlation function for the monomer positions. The essential physics in terms of the variational problem 
will be described by the Euler equations corresponding to free energy minimization. Although important 
conclusions concerning the overall size of the chain can be obtained |@] via the asymptotic analysis of the 
Euler equation, we will mainly pursue in this paper a numerical solution for finite number of monomers 
N. 

In section II we review the variational method as it was discussed previously for homopolymer chains 
and use a more generalized form to polyelectrolyte chain in theta and poor solvent |^-|^]. We will show 
how the different variational approach can be successfully applied to either open chains or chains with a 
cyclic geometry. The cyclic invariance || has the advantage to reduce the computational effort required 
to minimize our variational equations, because the chain coordinates can be decomposed in orthogonal 
(Rouse) modes. For open chains, the end effects have to be included explicitly. In section III present 
the results for cyclic chains of increasing size 2 4 < N < 2 16 . A global phase diagram in terms of solvent 
quality and inverse Bjerrum length is given. In section IV we consider the variational approach for open 
chains and we present the results for chains of length 2 4 < N < 2 9 . The open chain case result to be 
very involving, since the computational effort is enhanced of a factor N/ log 2 N with respect to the cyclic 
chain case, but inhomogeneous point to point correlation function are obtained and discussed. These will 
reveal the end effects we discussed for polyelectrolyte chains in theta solvent. Moreover, in the extended 
globule region, a non-monotonic correlation function, showing structural modulation along the chain was 
found, reminiscence of necklace formation. 




FIG. 1. For the two architectures considered, i.e., open chains a nd r ing chains, we expect for increasing values 
of the interaction strength, / 2 /b a first transition (see equation ( |f ,4| ) from globule to extended globule and a 
second one, at higher values of the interaction strength from extended globule to open chain (see equation (IFq). 



II. THE GAUSSIAN VARIATIONAL METHOD 

A. The homogeneous variational method - cyclic chains 

In this section, we apply the variational method (Gaussian approximation) to study the conformation 
of a polyelectrolyte chain in poor solvent. We compute the free energy of a self-interacting polymer 
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chain with Coulomb repulsion between the monomers, based on a discrete representation of a chain of N 
monomers. The method relies on the well known variational principle, where a Gaussian trial probability 
is proportional to the exponential of a quadratic form of the monomer coordinates fj . The variational free 
energy for a uniformly charged polymer of length TV in a d dimensional space can be computed according 
to the usual Gibbs Bogoliubov inequality 

F<F v = {H-H Q )o + F , (2.1) 



where Fo is the free energy of the Gaussian model defined by equation ( |A2|) . In equation (2.1), H 
represents the full Hamiltonian of a charged polymer chain with a long - range monomer-monomer 
interaction, 

v(f i -f-) = \f i -f j \- x , (2.2) 

and where is generally A < d. However throughout the present paper we use only the three dimensional 
Coulomb case, i.e., A = d — 2, d = 3, see also Appendix A. Following the Cloizeaux' || the expression for 
the Hamiltonian reads 

, N N N 

H/ksT = — 2 - nf - WY^m ri) + 

i—1 i—1 j—1 

N N N N N 

™/3! E E E s & - ^ft - + ^/ 2 /2 E E v ^ - ?i)> ( 2 - 3 ) 

i=l j=l k=l i=l j=l 

For the variational free energy Fy in the limit of long chain N ^> 1 , 

n—1 ^ ' n — 1 m— 1 

N , /i -[~7T 7 

+ u C h (3,l)^[6(n)]-i + — y log 5 (*)dfc--[l + log(27r/d)] ( 2 - 4 ) 



where b(n) = 1/b 2 ((fl — ri+ n ) 2 ) and is related to g(k) via equation ( A10| ) and (7^(3, 1) is a constant 



depending on the dimension d and interaction exponent A. We also express the coupling constants as 
dimensionless quantities u, r, w. The electrostatic interaction strength u is defined in term of Bjerrum 
length and charge fraction u = IrJ 2 /b. The second virial coefficient is replaced by t = v/b 3 . The 
third virial coefficient w is fixed to b 6 . During all numerical evaluations set w/b 6 = 1 for simplicity and 
normalization. Minimization of the variational free energy Fy yields 

g(k) = (1 - cob(*)) 
/ 3 \ i N 

- T (^j E( 1 - cos (" fc ))[ fe Hr l 

^ ' 71=1 

N N-n 

E {0- - cos(nk))[b(n)]~^[b(m)]-i + (1 - cos(mfe))[6(m)] _ 2[6(n)]"i} 

n—1 m— 1 

2iiM^(i_co8(n*))[6(n)]-i. (2.5) 



u 

6 „ 



Equation (2.5) is a self consistent equation for g(k). The variational principle we discussed is characterized 
by isotropic correlation functions, as we assumed (see Appendix A) considering Fourier representation. 
This is of course related to the cyclic chain architecture we assumed throughout the calculation. In the 
extended globular phase, onset of the Rayleigh instability, and the formation of necklaces, these isotropic 
and homogeneous correlation functions appear no longer sufficient to describe a structure, that, if shown 
to exist, is intrinsically inhomogeneous. Therefore we are going to generalize the above variational method 
for heterogeneous structures, releasing the translational invariance constraint. 
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B. Inhomogeneous variational method for an open chain 



In this section we generalize the variational method to study an open chain (without assuming cyclic 
invariance) and in d - dimensions. The free energy Fjsr depends on the entire spectrum of correlation 
function Q(i,j) = (flfj). Using the Gaussian trial Hamiltonian H Q = Q~ l {i, j)flfj we derive Fjy 

in theta solvent, 

, , , 2 N 

F N /k B T = --Tr log Q - -TV + — J] [Q(n, n) + Q(n + 1, n + 1) - 2Q(n, n+1)} 

n=l 

+ uC(d)^2 E [Q(n,n) + Q(m,m)-2Q(n,m)]-K (2.6) 

l<m <n<N 

where C(d) is a constant which depends only on dimensionality (see Appendix). In d = 3 for Coulomb 
interaction (A = 1), C(d) — (Sl 3 ) 2 tt /2) / (2n) 3 and SI3 = An is the solid angle in d = 3. The notation Q 
stands for matrix. Minimization of the free energy yields 

6Fn (2.7) 



SQ(i,j) 

and we obtain the Euler equation for an open chain corresponding to equation ( p.6[ ) in the form of a 
N x N matrix equation. 

d d^ 

^Q~\hj) = ^ p (hj) ~ ^C(d)J2[Q(i,i) + Q(n,n) - 2Q(i,n)]-iS ij 

U 3 

+ ^C(d)[Q(i,i) +Q(j,j) - 2Q(i,i)]-5(l - Sij), (2.8) 
P(iJ) = cn(i)5ij - S iJ+1 - S j:i+ i. (2.9) 

where matrix P(i,j) expresses connectivity in terms of the monomer coordinates r^. The ends of the 
chain satisfy different connectivity condition, cn(i) = 1 for monomers at the end of the chain i — 1, N 
and cn(i) = 2 for monomers along the chain. A general expression of the free energy for arbitrary 



solvent conditions is presented in Appendix B. The right hand side of Eq. 2.8 defines each element of 
the matrix Q -1 . This matrix equation has been solved self-consistently. The computational cost of 
matrix inversion scales like iV 3 . Note that the sum of each column (row) of Q^ 1 is vanishing (since the 
diagonal element Qu = — X^j^i Qij) indicating the linear dependence of the columns (rows). Due to the 
translational invariance the eigenvalue of matrix has a zero mode corresponding to the motion of the 
center of mass. This cause the determinant of the matrix to vanish. It is indeed necessary to introduce 
the constraint fjv = — 53< ^ to eliminate the degree of freedom corresponding to the motion of the center 



of mass. Proceeding to the first iteration of the self consistent Eq. 2.S, we define the reduced connectivity 
(N — 1) x (N — 1) matrix N r in the center of mass coordinate system. This can be obtained by expressing 
rV in terms of r^'s (i=l,..N-l), according to the following identity. 

N N N-1N-1 
i j i 3 

We obtain 

P r (i,j) = P(i,j) + P(N,N)-P(i,N)-P(N,j), (i = l,..N-l). (2.11) 
This operation still preserves the symmetry P r (i,j) = P r (j,i). We include the reduced form of the 



connectivity matrix P r (see Eq.2.8) and start the iterative scheme discussed above. Another way to 
remove the zero mode is by using the relative coordinate system (bond vector coordinates) bi = fi — fj_i 
which is the case described in f23|]. 



III. RESULTS - CYCLIC CHAIN 



In this section we present the results by solving the Euler equations (A13) and fl2.4|) for a polymer chain 



of length TV with cyclic boundary conditions. We compute the optimal profile g pt(k) that minimizes the 
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free energy Fy (see eq. (2.4)), for increasing values of the interaction strength u, in a theta solvent. At 
a given value of u we can obtain a direct measure of the blob size in the following way: we compute the 
real space correlation function b(n), related via Fourier transformation to the optimal profile g pt(k) (see 
equation ( Aldj )). In the weakly charged limit (i.e. for small values of u) the function b(n) — n 2l/ shows a 



crossover from Gaussian behavior (at small values of n, 
1) as described in equation (1 



v = 1/2) and a linear behavior (at larger values 
We define the blob size as the extrapolated crossover point 
of the exponent v between these two regimes. The result is shown in Fig. 2. 
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FIG. 2. Blob size, as measured from the direct solution of equations (A13) and (2.4), considering the crossover 



between the small distance behavior of the mean square distance b(n) versus the large distance behavior. Different 
values of the electrostatic blob size for increasing values of the interaction strength u are compared with the 
£ ~ l B behavior expected from phenomenological arguments. 



A. Neutral chain in Poor solvent. 



We now investigate the poor solvent regime, first considering a neutral chain. We will consider the 
effect of a Coulomb repulsion between the monomers in the next paragraph. In poor solvent conditions, 
as mentioned above, we expect the neutral chain to form a globule. There is one important length scale 
in the globule, namely the size of the density fluctuations. On length scales larger than the correlation 

1/2 

length £ T w bg T , the interaction between the monomers overcomes the thermal fluctuations, resulting 
in a dense, liquid like structure, made of Gaussian blobs. We were able to confirm quantitatively this 
scenario within our variational principle. We start from a Gaussian chain and let the solvent quality 
decrease. The density of the globule is determined by interplay of 2-body attraction and 3-body repulsion. 
As in the case of a theta solvent, the monomer-monomer correlation function shows a crossover between 
a Gaussian regime b(n) oc n and a regime where b(n) oc n 2 ' 3 . Moreover, to check that our variational 
method properly describes a neutral chain in poor solvent conditions we measure the overall size of the 
globule as a fu nction of the solvent quality. One should expect to recover results that are consistent with 
equation (Oj). 

In Fig. 3 the end to end distance is shown as a function of r. As expected the globule is formed as solvent 
quality becomes poorer (r increases) meanwhile the three body term is constant, as usually assumed in 
scaling theory. 
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Solvent quality x 



FIG. 3. End to end distance behavior in poor solvent for a cyclic chain of length N. The end to end distance 
is obtained for different r for chain length N=64, 128,256. The entire spectrum of correlation functions b(n) is 
measured at a given value of the solvent quality v ~ 0.4. 



In Fig. 4 the globule size is shown. The proper scaling behavior as a function of the chain length, 



predicted by scaling (see equation (1.3)), is recovered. 
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FIG. 4. Globular End to end distance R* versus chain length for chains of increasing length iV = 64, 128, 256. 
The continuous line is obtained by a fit of our numerical results, where the fit parameters are A — —0.9 and 
B = 0.3329 ± 0.0002. Similar results can be obtained for globules of different initial density. The proper scaling 
of R with density has also been checked carefully and confirms scaling results. 



We observe the overall size of the system to scale as R oc iV 1 / 3 . The same scaling behavior is ob- 
served when a weak enough Coulomb interaction is included within our variational approach. A sharp 
transition to extended globule appears for values of the Bjerrum length exceeding a critical value Zg as 
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shown in Fig. 5. According to the scaling theory of Khokhlov a sharp transition occurs due to Coulomb 
interactions between the monomers when the surface free energy becomes of the order of the free energy 
contribution due to excluded volume interactions. The variational method takes into account the surface 
terms responsible for the extended globule transition by a proper rescaling of the monomer-monomer 
correlation functions. 

Simple scaling arguments suggests that the critical value of the Bjerrum length where the transition 
occurs scales as 1/N where ./V is the number of monomers in the chain. We show in Fig. 5 the measured 
value of the critical Bjerrum length value as a function of the chain length N . Scaling properties of the 
end-to-end distance of the chain have been also discussed for the long - ranged intra- molecular interaction 
1 /r A at an arbitrary spatial dimensionality d . The validity of the variational approach for short range 
interacting homopolymer chains must be questioned J9| for natural reasons. In poor solvent the problem of 
logarithmic corrections that spuriously appear within the Gaussian Variational Principle cancels between 
the two and three body terms, so that a power law behavior for the end to end distance is measured, i.e. 
R oc TV 1 / 3 , as in Fig.4. 



B. Polyelectrolytes chains in poor solvent 



We now consider a polyelectrolyte in a poor solvent. When the polymer is charged the Coulomb 
repulsion between the monomers changes the shape of the globule but not its overall volume, defined only 
by the solvent quality. When the Coulomb repulsion -F cou i « kftTl-BpN 2 / R is comparable to the surface 
energy -F sur f ~ 7o^? 2 the globule either elongates into a cylindrical shape or splits into two connected 
globules. We expect an instability in the radius of gyration at finite value of u once we include a Coulomb 
repulsion in our vari ation al equations. Such a scenario is captured in the numerical implementation of 
the Euler equations (A1E). For a given chain length, the radius of gyration R has a sharp increase at a 
critical value of u*. According to Eq.1.4 we expect u* to scale as 1/N. This scaling behavior is shown in 
Fig.6. 
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FIG. 5. critical interaction strength u* versus chain length N. 



Fig.6 shows this instability for a given chain of lengths N — 64 and N = 128. The end to end distance 
and its derivative are computed as a function of u. 
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FIG. 6. End to End distance R for increasing values of u. The inset on the top left corner shows the derivative 
of R vs u. 



Following the sharp transition point in Fig. 6 allow us to construct the phase diagram for poor solvent 
conditions. The phase diagram, for different values of the interaction strength u, is shown in Fig. 7. 
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FIG. 7. Phase Diagram in poor solvent regime and inverse Bjerrum length, for a chain of length N = 2 . 
phase boundaries are represented by solid lines. The dashed line indicates that the transition becomes very weak 
according the hight of the second peak of Fig. 6. the two boundaries on the phase diagram discriminate: (a) the 
neutral globule phase, where the polyelectrolyte is under-charged and follows simple neutral chain statistics, (b) 
the cylindrical globule phase, where a multiple crossover for the monomer-monomer correlation function is ob- 



served from Gaussian (b(n) oc n, n < £ T ), to dense ( (b(r 



2/3 



N » n > £ T ) to extended (b(n) oc n ,n < N). 



(c) the stretched chain regime, where electrostatic blobs forms and the chain is fully stretched 
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In Fig. 7 we recover the different regimes described in the discussion of section I. The two boundaries on 
the phase diagram discriminate: (a) the neutral globule phase, where the polyelectrolyte is undercharged 
and follows simple neutral chain statistics, (b) the cylindrical global phase, where a multiple crossover 
for the monomer-monomer correlation function is observed from Gaussian (b(n) oc n,n < £ T ), to dense 
( (b(n) oc n 2 / 3 ,N >> n > £ T ) to extended (b(n) oc n 2 ,n < N). (c) the stretched chain regime, where 
electrostatic blobs forms and the chain is fully stretched (see Fig.2). 



IV. THE INHOMOGENEOUS VARIATIONAL METHOD 



A. Theta Solvent 



The application of Gaussian variational method to an open chain geometry allows us to investigate the 
end effect of finite chain which was neglected the cyclic geometry assumption of section III. In a theta 
solvent, the weakly charged polyelectrolyte has non uniform blob size distribution ]3l| indicating charge 
depletion at the center of the chain. In Fig. 8 we compute the size of the blob at different relative position 
of the chain. The estimation of blob size is obtained finding a crossover point from Gaussian to extended 
chain statistics as explained earlier. The dashed line shows the estimation from Eq.(1.2). 



logarithmic corrections in the form _/V(log-/V) 7 



1.4 



1.2 



o , 

o - ' 

o 



9' 



I I ■ ■ I ■ ■ i ' I 

-0.5 0.5 

Monomer position z/R n 

FIG. 8. Blob size £(z/R\\) at relative position z/R\\. The dashed line represents the theoretical prediction 
from Eq.1.2 with £m/R\\ = 0.25 



B. Poor solvent. Necklace formation 



In the cylindric globule phase, which we found (within the variational method for cyclic geometry) of 
section IIIB, we expect a correlation function that differs from the one computed in the globule phase. 
Any indication of globule splitting into sub-globules connected by strings will appear as a modulation 
along the chain of the monomer to monomer correlation function. The monomer-monomer correlation 
function b(i,j) = Q(i,i) + Q(J,j) — 2Q(i,j), in the globular phase, is almost constant in the middle of 
chain, as is shown as a lowest curve in Fig. 9. At the two ends of the chain, b(i,n) increases and this 
is an end effect in a similar way as discussed in the case of a chain in theta solvent. At higher values 
of the interaction strength u, we observe (see Fig. 9) two minima appears at both ends of the chain and 
the maximum at the center part of the chain. The minimum (maximum) in the correlation function 
corresponds to the density maximum (minimum). It is clear that for large value of it, the center part 
of chain is more stretched out and connects to two dense parts of the chain (pearls). The instability of 
the globule structure manifests clearly at the transition point u* ~ 0.24. The fluctuation of b(i,n) is 
relatively large as can be seen by the curve in Fig. 9. 
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50 100 150 200 

Monomer index i 

FIG. 9. The correlation function b(i,n + i), for various interaction strength u, is measured for fixed value of 
n = 7. The necklace formation proceeds the instability of the globular phase. 

The density maxima move to the edges of the chain with increasing u. We interpret this as a structure 
consisting of smaller pearls connected by strings. To investigate the detail structure of this new scenario, 
we plot in Fig. 10 b(i,i + n) for different value of n at the given interaction strength (u = 0.26 > u*) for 
a chain of length N = 200. For intermediate values of n, (n ~ 5), the density minimum clearly appears. 



3 

9 

15 



50 100 150 200 

Monomer index i 

FIG. 10. Monomer- monomer correlation function for values of the interaction strength u* < u < it** for a 
chain of length iV = 200. The figure shows b(i, n) vs i at different values of n = 3, 5, 7, 9, 15 



V. CONCLUSIONS 

We have shown in this paper that typical surface dominated physical effects in polymer globules can 
be derived beyond scaling theory. Variational techniques have been proven to provide satisfactory results 
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for systems with long ranged interaction potentials. Moreover, polymeric globules are systems where the 
density fluctuations inside the globules are small and are ideal candidates for variational computations. 
We have shown indeed a neutral globule can be described within our variational scheme, as well as 
the repulsive instability driven by Coulomb repulsion. The variational approach introduced in section 
IIA, and usually encountered in the literature to discuss a self interacting chain, suffers from the cyclic 
constrain that naturally hides end effects and structural modulation along the chain. For this reason 
we generalized the variational principle to the case of an open chain. The entire spectrum of monomer- 
to-monomer correlation functions was indeed obtained in this way. As mentioned above, the case of an 
open chain appeared to be more involved at the computational level, so that the chain length considered 
was one order of magnitude smaller compared to the case of a cyclic chain. Further improvement of our 
numerical techniques would eventually reveal a multiple splitting of the globule into sub-globules, beyond 
the simple "dumbell" formation we observed. 
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APPENDIX A: VARIATIONAL APPROACH UNDER CYCLIC INVARIANCE CONDITION 



Consider a Gaussian variational probability distribution 

Pv{n, ..,rjv) = Zy 1 exp{-H (ri, ..,r N )}, (Al) 

where Hq is given by 

, N N 

ff (rl, r N ) = - J2 E G (* - ? if ( A2 ) 

»=i j=i 

and Zy is a normalization constant determined by the condition 

P v (fi, ..,r N )d d n ■ ■■d d ? N = 1. (A3) 

By the appropriate choice of periodic boundary conditions for the chain, the quadratic form Hq can 
be diagonalized by introducing cyclic coordinates (Rouse modes) for the monomer positions 

JV 



p q = N 2 exp[ i2njq/N]r : j 
3=1 

N 

r-j = N~i ^2 ex Ph i27T jq/N}p q . (A4) 

Due to cyclic invariance the Cartesian components of p q satisfies 

MPpT) = ^f5 qq ,g-\2ir q /N), (A5) 

where the function g^ ! L (2nq/N) is posi tive and discontinuous function of k = 2nq/N, and is related to 
the Gaussian propagator G of equation (A2) via the following expression 

N-l 

g(2nq/N) = ^ G{n)[l - cos(2vrgn/A r )], (A6) 

71=1 

The variational free energy Fy in the limit of infinitely long chain N — > oo, 
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71=1 



d_ 

-in 



logg(k)dk--[l+log(27r/d)], 



where 



F h [g(k)} 

k B TN 



[b{r 



(A7) 



(A8) 



n—l m—1 



and Ch(d, A) = C(d, A) x d? . For arbitrary A, 



C(d, A) = 2' _1 T(d - A - l)r(-) sin 



n(d- A-l) 



(r» d ) 2 /(2^) d 



where fid = 27T2 /T(4) is the solid angle for d-dimensions. For Coulomb interaction A = d 
simplifies to C(d) = ^/(2vr) 2 ^2. 



The function g(k) is related to b(n) by 



being, 



The function <?(fc) satisfies the symmetry requirements 

g{k)=g{-k), g(k + 2n)=g(k). 
Minimization of the variational free energy Fy yields 



(A9) 
2,C(d) 

(A10) 

(All) 
(A12) 



;(Jb) = (1 - cos(fc)) + u-C h (d) V(l - cos(nk))[b(r 

ft < J 



— T 



<f\* N 



2n 



— cos(nfc))[&(r 



71=1 
d iV iV— 771 



+ u- — £ ^ {(l-cos(nfe))[6(n)]-^[6(m)]-5 +(i-cos(mfe))[6(m)]-^ 2 [6(n)]-^ (A13) 



' n=l 771=1 

which reduces to equation (E75J) for A = 1 and d = 3. 



APPENDIX B: THE VARIATIONAL METHOD FOR THE INHOMOGENEOUS OPEN 

CHAIN 



In this appendix we write explicitly the Euler equations for the inhomogeneous variational method, 
including terms corresponding to two and three body interactions. In this case the monomer-to-monomer 
correlation function Q(i,j) depends from both its arguments, not just from the distance i — j, as in the 
cyclic chain case. It is convenient to choose the trial Hamiltonian Ho in the form 

1 N N 

H (f 1: . . . ,f N ) = - £ £ Q~ X U, iWi (bi) 

3 = 1 1=1 

The straightforward calculations leads to the following result for the variational free energy: 
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Fv(Q) = -^TrlogQ --N+—J2 b(n, n + 1) + -C(d) E E [&( 



n=l m n^m 

- ■' 

T I 1 \ 2 X - ^ X - ^ r, / m-^ W f 1 



2 V2tt 



EE t^™)r 



6 V27T 



EE E [K*,n)] _I [&(*."»)]■ 

k n^k m^k^n 



where 



b(n, m) = Q(n, n) + Q(m, m) — 2Q(n, m) , 



(B2) 



(B3) 



As in the homogeneous variational method, we proceed to compute the Euler equations, obtained mini- 
mizing the above expression for the free energy in terms of the propagator Q(i,j). 



7—zF v {Q(i,j)} = Q 



Given that 5 
d 



TrlogQ 



5Q(i,j) 

/SQ(i,j) = Q~ 1 (i,j) one gets 



Td / 1 



u><i / 1 



AT 



rd / 1 



+ Y i) E[K^)]--^-t(^) [^i)]--d-^) 



2 V 2tt 



+ 



6 V27T 



u>c? / 1 



E E '[K*."»)] 



/ 1 



6 V2tt 



d N 



6 V27T 



(B4) 



^[KiJT^ {[b(i,m))]-t + [bV,m))]-*}(l-5 ij ). 



(B5) 
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